ASYMPTOTIC AND BIFURCATION ANALYSIS OF WAVE-PINNING IN A 
REACTION-DIFFUSION MODEL FOR CELL POLARIZATION 
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Abstract. We describe and analyze a bistable reaction-diffusion (RD) model for two interconverting chemical species 
that exhibits a phenomenon of wave-pinning: a wave of activation of one of the species is initiated at one end of the 
domain, moves into the domain, decelerates, and eventually stops inside the domain, forming a stationary front. The 
second ("inactive") species is depleted in this process. This behavior arises in a model for chemical polarization of a cell 
by Rho GTPases in response to stimulation. The initially spatially homogeneous concentration profile (representative of a 
resting cell) develops into an asymmetric stationary front profile (typical of a polarized cell). Wave-pinning here is based 
on three properties: (1) mass conservation in a finite domain, (2) nonlinear reaction kinetics allowing for multiple stable 
steady states, and (3) a sufficiently large difference in diffusion of the two species. Using matched asymptotic analysis, we 
explain the mathematical basis of wave-pinning, and predict the speed and pinned position of the wave. An analysis of 
the bifurcation of the pinned front solution reveals how the wave-pinning regime depends on parameters such as rates of 
diffusion and total mass of the species. We describe two ways in which the pinned solution can be lost depending on the 
details of the reaction kinetics: a saddle-node or a pitchfork bifurcation. 
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1. Introduction. In a recent reaction-diffusion (RD) model for biochemical cell polarization pro- 
posed in [21] we found a wave-based phenomenon whereby a traveling wave is initiated at one end of 
a finite, homogeneous ID domain, moves across the domain, but stalls before arriving at the opposite 
end. We refer to this behavior as wave-pinning. We observed that this phenomenon was obtained from 
a two-component RD system obeying a modest set of assumptions: (1) Mass is conserved and limited, 
i.e. there is no production or removal, only exchange between one species and the other. (2) One species 
is far more mobile than the other, e.g. due to binding to immobile structures, or embedding in a lipid 
membrane. (3) There is feedback (autocatalysis) from one form to further conversion to that form. 

The biological motivation for studying our specific system comes from internal chemical reorganization 
that is the initial stage of polarization of a living eukaryotic cell, such as a white-blood cell, amoeba, or 
yeast in response to a signal. Such chemical asymmetry then organizes the downstream response of the 
cell (e.g. shape change, motility, division, etc). Explaining the basis for such symmetry breaking has 
become an important question in cell biology over the past decade, motivating such mathematical models 
as [201 [Ml HS1 [5] ■ Our own work has focused on the role of switch- like polarity proteins called Rho 
GTPases, which are conserved in eukaryotic cells from amoebae to humans. Upon stimulation, levels 
of Rho GTPase activity rapidly redistribute across a cell. For example, some members of this family 
(Rac, Cdc42) become strongly activated at one end (which subsequently becomes the front of the cell 
[16j [23]) whereas others (such as Rho A) dominate at the opposite end (which becomes the rear [43]). 
Whereas in our previous work we investigated such phenomena in the context of the actin cytoskeleton 
and cell motion, |18[ [5] , here we are concerned only with the mathematical basis for the initial symmetry 
breaking. Originally, we explored multiple interacting Rho GTPases, to determine how interactions 
between several members of this family affect spatio-temporal dynamics (181 112] . In the more recent work 
|21j . we investigated a minimal system, consisting of a single active-inactive pair of GTPases. From a 
mathematical perspective, this yields an opportunity for deeper analysis. From a biological perspective, 
it clarifies what are minimal conditions required for symmetry breaking. 

The model described here and in [5T] is consequently based on the following abstraction of exper- 
imental observations about Rho proteins: (1) The protein has an active (GTP-bound) and an inactive 
(GDP-bound) form. (2) The active forms are exclusively found on the cell membrane; those in the fluid 
interior of the cell (cytosol) are inactive. (3) There is a 100-fold difference between rates of diffusion 
of cytosolic vs membrane bound proteins [55]. (4) Continual active- inactive exchange is essential for 
proper polarization. If this exchange is stopped, the cell cannot polarize [TU]. (5) On the time-scale of 
polarization (minutes), there is little or no protein synthesis in the cell (timescalc of hours), i.e. during 
polarization, the total amount of the given protein is roughly constant. (6) Feedback from an active form 
to further activation are common. A schematic diagram of our model is given in Fig. 12.11 Cases where 
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Fig. 2.1. Left: The cell (initially unpolarized) consists of membrane (shaded) and cytosol (white). The ID model 
112.11 1 represents the chemical distribution of two forms of a protein (u, and v) along the diameter of a cell, idealized as a 
thin flat strip of uniform thickness. (The nucleus and other nonuniform features are neglected). Right: in the RD system, 
u(x,t) is an active protein, resident in and diffusing slowly along the membrane, v(x,t) is an inactive protein diffusing in 
the cytosol. Interconversion is subject to positive feedback from u to itself. The cell diameter is along < x < L. Cell and 
compartment sizes not drawn to scale. 



this has been established experimentally include HOI E3 HZ] ■ 

The RD system in this study encapsulates all the above aspects, and exhibits self-polarization as 
a result of the wave-like phenomenon described above. The purpose of this paper is to investigate the 
properties of this model. In Section [2j we formulate the model in one space dimension. We first seek to 
obtain a mathematically clear picture of wave-pinning. This is achieved by way of matched asymptotic 
calculations, presented in Section [3] Here, we shall see how the wave speed, shape and stall positions are 
affected by the parameters of the problem. We also briefly discuss the behavior of higher dimensional 
generalizations of the one-dimensional model. Next, we study the bifurcation structure of our system in 
[4] We delineate the parameter regime for which wave-pinning is possible, and describe the bifurcation 
structures that are possible for different reaction kinetics. A summary and biological implications are 
presented in the Discussion. 

2. Model formulation. Consider a one dimensional domain il = {x : < x < L}. Denote by 
u{x,t) and v(x,t) the concentrations of active and inactive protein respectively at position x and time t. 
This one-dimensional model would be valid for a flat cell of sufficiently small thickness so that appreciable 
chemical gradients do not develop in the thickness direction (Fig. I2.1[) . In this case, both membrane and 
cytosolic positions can be described by a single coordinate x, and thus, we may treat the membrane-bound 
species u and the cytosolic species v as residing in the same domain f2. There are biological situations 
that warrant models in higher spatial dimensions, and we will consider such generalizations in Section 
13.41 From a mathematical point of view, however, we shall see that much of the behavior of interest is 
already present in the one-dimensional model. 

The concentrations u and v satisfy the following equations 



du d 2 u 



Of 

dv 



' dx 2 
d 2 v 



dt= D »dx^- fiu > v) > 



(2.1a) 
(2.1b) 



where f(u, v) is the rate of interconversion of v to u, and the rates of diffusion satisfy D u <C D v , reflecting 
the fact that the membrane bound species u diffuses much more slowly than the cytosolic species v. The 
boundary conditions are 



du dv 



dx dx 

It is clear that system (|2.1[) leads to mass conservation, 



= 0, x = 0,L. 



(U + v)dx = A'total, 



where Atotai is a time- independent constant. 

The following reaction term f(u,v) was proposed in [21j : 



(2.1c) 



(2.2) 



f(u,v) = T) [5 + 



v — r/u 



(2.3) 
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where 77, 7, m > 0, S > are constants. The above reaction term is written as the difference between a 
production and a decay term. The production term can be seen as v times the production rate. The 
production rate has a sigmoidai shape as a function of it, which expresses the presence of positive feedback 
[33JI3H]- For suitable choices of 7, m and S, f(u, v) has the following property. The expression f(u, v) = 0, 
seen as an equation for u with v fixed over a suitable range, has three roots U-(v) < u m (v) < u+(v). 
Moreover, u±(v) arc stable fixed points of the ODE 4| = f(u,v) whereas u m (v) is an unstable fixed 
point. In other words, the function f(u,v) is a bistable function of u over a range of v values. Much 
of the analysis to follow applies not only to the specific form of f(u,v) given in (|2.3p but to a family of 
reaction terms satisfying a number of properties including bistability. A precise characterization of this 
family will be given shortly. 

We now make our equations dimcnsionlcss. We scale concentrations with m and the reaction rate 
with 77, both of which are dictated by the form of the reaction term (see (|2.3[) ). Take the domain length 
L to be the relevant length scale. Equations (|2.ip can be rescaled using 

L , , 

u = mu, v = mv, x = Lx, t = t, (2-4) 

VV D u 

where u, v, x, and t are dimensionless variables. The scaling in time is chosen so that we obtain a 
distinguished limit appropriate for the analysis of wave-pinning (see next Section). We define: 



,2 



e 



Da n D v 



D = -T2- ( 2 -5) 



rjL 2 ' 77L 2 



Given D u <C D v , we let e be a small quantity. Wc let D = 0(1) with respect to e. This assumption may 
be written as \J D v /r/ ss L, i.e. on the timescale of the biochemical reaction, the inactive substance can 
diffuse across the domain. In the context of cell polarization, we have a typical cell diameter L ps 10/im, 
reaction timescale r\ rj Is -1 , and diffusion coefficients D u = 0.1 /.im 2 s _1 and D v = 10/.im 2 s _1 . The 
dimensionless constants are then e sa 0.03 and D as 0.1. 

Substituting the relationships (|2.4j) and (|2.5[) into (|2.ip dropping the ~ and using the same symbol / 
for the dimensionless reaction term, we obtain: 

£ S = e V + fM > (2 ' 6a) 

dv d^v 

with boundary conditions: 

|H = g = „, .. M . ( ,ec, 

Note that our domain is now < x < 1. The reaction term (|2.3I) assumes the following dimensionless 
form: 

f(u,v)= (s+^^jv-u. (2.7) 

The (dimensionless) total amount of protein satisfies 

1 

(u + v)dx = K. (2.8) 

where K = K totll \/m. Wc shall henceforth work almost exclusively with the dimensionless system. 

As mentioned earlier, we shall consider not only (|2.7p but a family of reaction terms satisfying the 
following properties: 

1. (Bistability Condition) In some range v m - m < v < u max (bistable range), the equation f(u, v) = 
has three roots, U—(v) < u m (v) < u + (v). Keeping v fixed within the bistable range, u±(v) are 
stable fixed points and u m (v) is an unstable fixed point of the ODE ^ = f(u,v). That is: 

df df 

^(u±(«),«) < 0, ^(u m (v),v) > 0. (2.9) 
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2. (Homogeneous Stability Condition) The homogeneous states, (it, v) = (u±(v), v), v m - m < v < v max 
are stable states of the system (|2.6|) . 

3. (Velocity Sign Condition) There is one value v = i> c ,u m i n < v c < v max at which the following 
integral I{v) vanishes: 

/""+(") 

I(v) = / f(u,v)du. (2.10) 

J U- (v) 

We assume in addition that / > for v > v c and / < for v < v c . 
The first condition is the bistability condition that was mentioned earlier. The reason for the name of 
the third condition will become clear in the next Section. We shall see in Section 13.11 that the second 
condition can be reduced to the following: 

<0. (2.11) 

(u,v) — (u±(v),v) 

Assuming this result, we can check that (|2.7p satisfies the above properties for the following parameter 
values. For 7 > and 6 > 0, (|2.7p satisfies the above conditions if and only if 

7 > 86. (2.12) 

The corresponding bistable range is given by v m - m = k+ < v < «_ = u max where: 









dv J 



K± = ;U + TT^J ' W± = V 2(1 + p) = ? (2 - 13) 

When 5 = 0, v m \ n = 2/7 and v max — 00. In our computational examples, we shall make use of (|2.7p with 
5 = and 7 = 1, which we record here for future reference: 

/(«>«) = I^-«- ( 2 - 14 ) 
In this case, uq(v) and u±(v) can be computed explicitly: 



V — V — 4 V + "yV^ — 4 
u-(v) = 0, u m (v) = ,u+(v) = . (2.15) 

We shall often make use of the following caricature of (|2.7|l : 

f(u,v) = u(l-u)(u-l-v). (2.16) 

It is easy to check that (|2.16p satisfies all of the above properties. For this reaction term, the bistable 
range is < v < 00. This example makes certain algebraic manipulations easier than (|2.7p or (|2.14p . In 
Section we shall also make use of another cubic that satisfies the above conditions: 

f{u,v) = -(u- l)(u-u m )(u + 1), u m = a > 0. (2-17) 

yjl + {av) 2 

The bistable range for the model with kinetics (|2.17p is — 00 < v < 00. At least one of the roots 
of this polynomial is always negative, and thus, it is no longer possible to interpret u and v as being 
concentrations of chemicals. The arguments to follow, however, never require that u and v be positive. 
Both (|2.16p and (|2.17p will prove useful in understanding the bifurcation structure of our system. 

We now describe the behavior that we wish to explain. If we consider (|2.6a[) as a stand-alone equation 
for fixed v, it is a scalar reaction diffusion equation of bistable type. It is well-known that such equations 
support propagating front solutions when posed on an infinite domain. Coupling this with (|2.6b[) on a 
finite domain gives rise to wave-pinning. In Fig. 12.21 we show simulation results for our dimensionless 
system (|2.6p with the reaction terms (|2.14p and (|2.16p . The concentrations are initialized so that u is 
high close to x = whereas v is spatially uniform. This represents a stimulus at the left end of the 
domain. The initial rectangular profile of u develops into a steep front which propagates into the domain. 
The height of the front gradually changes and the front eventually comes to a halt. The left portion 
of the domain has a high concentration of the active species u whereas the right portion of the domain 
has a low concentration. The spatially localized initial stimulus has been amplified to produce a stable 
spatial segregation of the domain into a "front" and a "back". The one-dimensional cell has achieved 
polarization. 
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Fig. 2.2. Wave-pinning behavior for the reaction diffusion model (1 2 . 6 1> with parameters e = 0.05, D = 1. (a) Hill 
function reaction kinetics (12. 14ft with S = 0, -y = 1, m = 1, K = 2.8. (b) Cubic reaction kinetics 112. 1611 and K = 1.9. 
Solutions to u (solid) and v (dashed) are shown at the indicated times. The wave is initiated as the square pulse in u, as 
shown at t = 0. 



3. Asymptotic Analysis of Wave-Pinning. In this Section, we perform an asymptotic analysis 
of wave-pinning to obtain the speed of the wave and its stall position. We shall first deal with the ID 
RD system (|2.6[) and consider its higher dimensional generalizations in Section [ 



3.1. Stability of the Homogeneous State. Let (u s (v),v) be a steady state of (|2.6p . where 
u s = u± or u m . Linearize (|2.6p about (u s ,v): 



^( u ]-c( u ]=j( u ]+^( e2u ) J=( f » 



(3.1) 

(u,v) — (u s (v) ,v) 



where f u and /„ denote partial derivatives of / with respect to u and v respectively. Here, the Jacobian 
of the reaction terms, J, is evaluated at (u s , v). To study linear stability, we study the spectral properties 
of the operator C under boundary conditions (|2.6c|) . We must also respect the mass constraint (|2.8|) so 
that the perturbations satisfy: 

1 

(u + v)dx = 0. (3.2) 

Since we are on a bounded domain, we need only consider eigenvalues. We thus consider the eigenvalue 
problem: 

'(:)-*(:)■ m 

where u,v must satisfy (|2.6cp as well as (|3.2j) . It is clear that all eigenfunctions are of the form: 

coskx (3.4) 



v [a. 



where k = rvn, n = 0, 1, 2, • • • where a u and a v are constants such that {a u , a v ) ^ (0, 0). When k ^ 0, 
a u and a v are arbitrary, whereas when k = 0, a u +a„ = to satisfy (|3.2[) . Is is easily seen that the 
eigenvalues satisfy the quadratic equation: 

A 2 - r k \ + A k = 0, r k = tr£ k , A k = detC k , C k = y^ f * fu _ D tf_ ;J ( 3 - 5 ) 

where tr£^ and det£fe denote, respectively, the trace and determinant of the 2x2 matrix C k . Let us 
first consider the case k = 0. In this case, the two solutions to the above quadratic equation are: 

A = 0, or A = t = f u - /«. (3.6) 

If To < 0, then the second eigenvalue is negative. As the reader can easily check, as soon as we assume 
t = f u — f v ^ 0, the cigenfunction associated with A = ceases to satisfy the mass constraint. Thus, if 
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we assume tq < we have stability for k = 0. Now we turn to the case k ^ 0. Both roots of f|3 . 5[) have 
negative real part if and only if T& < and Afc > 0. Since 



Tfc < so long as t < 0. 



r fc = -(D + e A )k 2 +T <r (3.7) 



A fe = £e 2 fc 4 - f u (D - e 2 )k 2 - r Q e 2 k 2 . (3.8) 



Therefore, Afc > so long as /„ < and D > e 2 . The condition D > e 2 is always met since we are 
assuming that e is small. For u s = u±, f u < is met by the bistability condition (|2.9[) . Thus, for / 
satisfying the bistability condition, the homogeneous stability condition of the last Section is equivalent 
to To < 0. This is condition (|2.1ip . It is interesting that the stability condition of the ODE system with 
fixed v (f u < 0) together with the stability condition for spatially homogeneous perturbations (to < 0) 
implies stability for all wave numbers. 

For u s = u m , f u > by (|2.9[) . For fixed k, (|3.8[) can be made negative by making e sufficiently small, 
and thus (u m (v),v) is always an unstable steady state for small enough e. This does not preclude the 
possibility that (u m (v),v) be a stable steady state for some finite e value. Suppose f u > and tq < 0. 
Let us consider the positivity of Afc, k > ir: 

|| = De 2 k 2 - f u (D - e 2 ) - r e 2 > De 2 ir 2 - f u (D - e 2 ) - r e 2 . (3.9) 

Therefore, (u m (v),v) is a stable steady state of the system so long as the right-most quantity is positive. 
This is the case if e satisfies the following bound: 

>' > d£tt„- < 3 - 10 > 

Since To < by assumption, /„ < f v and thus the right hand side of the above inequality is less than 
D. Therefore, if To < 0, there is a range of values satisfying e 2 < D (i.e. diffusion coefficient of u is 
smaller than that of v) for which (u m (v), v) is a stable steady state of (|2.6I) . On the other hand, if To > 0, 
(u m (v),v) is always unstable. 

For (|2.16[) , it is easily seen that (u m (v),v) is always unstable since /„ = at u = u m = 1 and thus 
r o = fu > 0. For (|2.17p . (u m (v), v) can be stable for a range of e values if a > 1 and v is in a suitable 
range. The stability of this middle stationary state does not play a role in the wave-pinning analysis of 
the next Section. However, it does play a role in determining the bifurcation structure of the system as 
we shall see in Section [Ol 



3.2. Asymptotic Analysis of Wave-Pinning. We now consider the dynamics of (|2.6p . There 
are three time scales in this model, the short, intermediate and long time scales. The intermediate time 
scale is of greatest interest to us, and (|2.6[) is scaled accordingly. We shall start with a brief discussion of 
the short time scale. Discussion of the long time scale will be deferred to the next Section. 

We introduce the short time variable t s = t/e. Then, (|2.6[) may be rewritten as: 



du 



dt* dx 2 



+ f(u,v), (3.11a) 



dv d 2 v 

m s = D ^- fM - (3 - ub) 

with no-flux boundary conditions. Assuming that u admits an expansion in e of the form u = uq + eu\ 
(and likewise for v), and substituting this into the above, we find that uq and i>o satisfy the equations: 

^- = f(u ,v ), (3.12a) 
ot s 

Suppose t>o satisfies w m i n < Vq < v max so that f(uo, vq) is bistable in uq. The first equation tells us that uq 
will evolve towards either u + (vq) or u_(«o) depending on whether uq(vq) is greater or less than u m (vo)- 
At the end of the short time scale, Vq will have a spatial profile that is uniform whereas uq will assume the 
values of u + (vq) or u_(t>o) depending on position. In other words, the domain will have segregated into 
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regions where uq = u+(«o) or w_(i>o). This profile will serve as our initial condition for the intermediate 
time scale. In general, this initial profile consists of multiple transition layers where u switches its value 
from u+ to it_ or vice versa. If there are no transition layers, this is nothing other than the stable steady 
state whose stability we just studied. In this Section, we shall restrict our attention to the case when the 
initial profile consists only of a single transition layer. Analysis in the case of multiple transition layers 
is essentially the same, and will be discussed in the next Section. 

We now begin the analysis in the intermediate time scale. Let <j>(i) be the position of the transition 
layer or the front. Note that the the position of the front changes with time. We now perform a matched 
asymptotic calculation. 

Expand u = uq + eu\ ■ ■ ■ and likewise for v. Substituting these expansions into (|2.6b .b) and retaining 
leading order terms we have the following equations for uq and vq: 

= f(uo,v o ), (3.13a) 

= D^-f(u Q ,v o ). (3.13b) 

Equations (|3.13[) are valid in the outer region < x < 4>{t) — 0(e) and <fr(t) + 0(e) < x < 1, that is, 
at some small distance away from the sharp transition zone at the front. Note that it is impossible to 
solve the above system with most initial data for uq and Vq. This is the reason why we need to insert 
a short time scale before this intermediate time scale. During the short time scale, the arbitrary initial 
condition evolves into an initial profile that is admissible as an initial condition for the intermediate time 
scale analysis. Adding p.lSb .b 1 ). we find: 

From (|3.14[) and boundary conditions (|2.6c[) . we conclude that: 

i*(x, t ) = H*> **«m- 0(e), 

1 ' \v>(t) <f>(t) + 0(e)<x<l, 1 ' 

where the values of v to the right and to the left of the front, v > and v<_, do not depend on x. From 
(|3.13aj) . uq takes on one of the values u+,u_ or u m in the outer regions. Since we are seeking a front 
solution, we let: 

u (x,t) = H;<l <?(*)-<>(*), (3.16) 

0V ' \u_(»>) <f>(t) + 0(e) < x < 1. V ' 

We have assumed, without loss of generality, that u transitions from u + to «_ in the direction of increasing 
x as we traverse <f>(t). 

Let w(x,t) — x — (j)(t) be the distance from the transition front, i.e. w = at the front. Introduce a 
stretched coordinate £ for the inner layer close to the evolving front: 

e= ^ = ^iW. (3.17) 

The inner solution is denoted by U, V, where 

U(£,t)=u((x-<j>(t))/e,t), V(£,t)=v((x-4>(t))/e,t). (3.18) 

Note that (|3 . 1 8[) is not a traveling front solution in the strict sense, as the wave speed d<p/dt is not 
constant. As the amplitudes of U and V also change with time, we do not assume u{x,t) = E/(£), but 
rather u(x, t) = U(£,t), and likewise for V. 

Substitute the new scaling into (|2.6[) to obtain the inner equations on — oo < £ < oo: 

f(U,V), (3.19a) 
-f(U,V). (3.19b) 



dU 


d<p dU 


d 2 u 




~di~dti ~ 


&e + 


dV 


d<j> dV 


Dd 2 V 




dt d£ ~ 


e 2 
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Expanding U, V and <p in powers of e, we obtain, to leading order, 

d 2 U d(j)Q dU 



d£ 2 dt 9£ 

d 2 v 



+ /(t/ ,Vb)=0, (3.20a) 



From (|3.20b|) . it follows that 



0. (3.20b) 



V = ai(f)£ + a 2 (i), (3.21) 



where ai(t), 02 (i) are arbitrary functions of t to be determined from matching. 
We match the inner (Vo) and outer (vo) solutions, 

lim V (£)=v < , lim Vo (£) = «>• ( 3 - 22 ) 

For these limits to exist, Vo must be a constant in the inner layer, i.e. 

v = V . (3.23) 

Thus, Vo is spatially uniform throughout the domain, and is equal to the outer solution vq. We thus 
recover our observation that vq should be uniform by the time the dynamics in the intermediate time 
scale dominates. We drop the dependence of vo on x (and Vq on £). 

We next consider a solution for Uq in the inner layer. Since Vo is spatially constant in the inner 
layer, (|3.20ap is an equation in Uq only, where Vo is a parameter (that varies in time). We must solve 
the boundary value problem (|3.20aj) with the matching conditions from (|3.16p as boundary conditions at 
±oo: 

lim U (0 = u+(V ), lim U (0 = U-(V ). (3.24) 

£— 00 £— too 

Such a heteroclinic solution Uq(£, Vq), unique up to translation, exists for general bistable reaction terms 
f(U, V) [HH1]. Multiplying ()3.20a|) by dU$/d£ and integrating from £ = -00 to £ = 00, we obtain: 



r + %°h(s,Vo)ds 
c(Vo)= y Vo) \ 2 ■ (3.25) 



dt 



An explicit analytical expression for c(y) cannot in general be obtained. An exception is when the reaction 
kinetics is of the form /(it, v) = —(u — u + (v))(u — u m (v))(u — it_(u)), where u_ < u m < u+. In this case 
c(v) is given by }2"2"] : 

c(v) = -j= («+(«) - 2u m («) + u_(»)) . (3.26) 

The sign of the velocity, however, is determined by the numerator of fraction in (|3.25p and can thus be 
easily determined given the reaction term f(u, v). By velocity sign condition (see equation (|2.10p ) we see 
that d<po/dt is positive when Vq > v c and negative when Vq < v c . 
By (|2.2p . we see that uo and Vq = Vq satisfy the relation: 

vq + / uodx = K. (3.27) 
Jo 

The integral of uq can be approximated by contributions from the two outer regions (to left and right of 
the front) and a 0(e) contribution from the inner region: 

1 M)-<D{e) ,1 

Uodx= uodx+ Uodx + O(e) 

JO J<l>(t)+0(e) 

= u+(v )Mt) + u_(«o)(l - ^o(*)) + O(e), 

where we have used (|3.16p in the second equality. Discard terms of 0(e). The reaction-diffusion system 
is then reduced to the following ordinary-differcntial-algcbraic system: 

-J^=c(w ), Vq = K -u + (v )(j)o - U-(v )(l - fo), (3.28) 
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where c(v ) is given by (|3.25f) . In (|3.28[) . the total amount of material, K, is allocated to a band of width 
cf>o at level u+, a band of width 1 — </>o at level u_ , and a homogeneous level of vq across the entire interval. 

We now show that the front speed, defio/dt, and the rate of change dvo/dt, have opposite signs. 
Differentiating the relation f(u±(v),v) = with respect to v and using (|2.1ip leads to 



= 

Using (|2.9[) we conclude that: 



df du ± df_ 

du dv dv 



J.±(v) 



du±\ df 
dv J du 



(3.29) 



du± > 0. (3.30) 



dv 

Differentiating the second relation in (|3.28[) with respect to t results in: 

+ — -j 0o + — -, (1 -<h))-jr = -( u +( v o) - u-(v )) — . (3.31) 

dv dv J dt dt 

Since the front position must reside within a domain of unit length, we have < </>o < 1. Using this and 
(|3.30[) . we see that the factor multiplying dvo/dt in (|3.31[) is positive. Since (u + — U-) > 0, we conclude 
from (|3.3ip that dvo/dt and d<j)/dt have opposite signs. Thus, vq is depleted as the wave progresses across 
the domain. It is interesting that this conclusion was obtained using the two conditions, bistability and 
homogeneous stability. 

Suppose v is sufficiently large initially, i.e., vq > v c at t = 0. Since d<f>o/dt is positive for Vq > v c , 
dvo/dt < 0. Thus, vq decreases as the front </>o advances. If vo approaches v c the front will come to a 
halt, i.e. will become pinned. Suppose the front is pinned at <f> p . Then </> p can be determined as follows. 
When the wave pins, we have 

v c = K - u + (v c )4> p - ti_(u c )(l - 4> p ). (3.32) 

We can interpret p.32[) as a relation between <j> v and K. We must have < <j) p < 1. This leads to a 
condition on K for wave-pinning to occur: 

v c + m_(u c ) < K < v c + u + (v c ) (3.33) 

that is, for wave-pinning to occur, the total concentration of chemical in the domain must fall within a 
range given by (|3.33[) . The pinned front is stable; if the front is perturbed, it will relax back to the pinned 
position <p p as can be seen from the velocity sign condition and the fact that ^f- and ^ have opposite 
sign. 

We now illustrate the above theory with the reaction term (|2.16[) . In this case, the reaction term is 
a cubic polynomial in u, and we may apply (|3.26p to find an explicit expression for c(v). The leading 
order equations (|3.28[) become: 

d(j> Q vp-1 . . 

— = — -j=-, v =K-{l + v )(j}o. (3.34) 

From (|3.34p . we find that the wave stops when Vq = 1 = v c . Condition (|3.33p reduces to: 

1 < K < 3. (3.35) 

Solving (|3.34p for vo, we obtain 

*p.= l(£4-l), v = ^l. (3.36) 
dt V2 V 1 + / 1 + <f>o 

The position at which the wave stalls, is therefore 

0p = (3.37) 

Fig. 13.11 shows that predictions of the ODE (|3.36p agree with numerical solutions to the full PDE system 
(|2.6p using the cubic reaction kinetics, (|2.16p . The exact front position is calculated from the numerical 
solution of the PDE system by tracking the position num at which u = u m (v) (u m = 1 for reaction 
kinetics (|2.16p ). n um(io) is used as an initial condition, where to ~ is a time at which the solution to 
the PDE system has relaxed to the form assumed in the asymptotic calculations. The error decreases 
with time as the wave becomes pinned. Based on the numerical evidence, we find that the leading order 
approximation is accurate to order e. To get a measure of the error of the leading term approximation, 
we can calculate the next term in the asymptotic expression. We refer the reader to [TT]. 
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(a) (b) (c) 

Fig. 3.1. (a) The evolution of the front position (f> n um obtained by numerically solving the full PDE system (12.66 with 
reaction kinetics 112.161 1 (solid), and the solution of the zero order asymptotic order approximation, (f>0 from Eq. A3. 361 1 
(dashed), (b) The CTTOT (pnum — 4>0 over time (c) The effect of e on the &TTOT (pnum — 

00. Parameters: e = 0.05, D — 1. 



3.3. Multiple Layers and Long Time Behavior. In the previous section, we discussed the 
behavior of system (|2.6[) in the intermediate time scale under the assumption that the initial profile 
consists only of a single front. We discuss what happens when the initial profile has multiple fronts or 
layers. Let 4>k(t), k = 1, ■ • • n be the front positions so that 4>k(t) < (f>k+i(t). For notational convenience, 
we let <po(t) = and 4> n +i{t) = 1. If u transitions from w + to u_ as we cross a front in the positive x 
direction, we shall call this a positive front. If the transition is from u_ to u + , we call this a negative 
front. In the sequel, we shall assume that 4>\{t) is a positive front. The case in which 4>i(t) is a negative 
front can be treated in an analogous fashion. If (f>i(t) is a positive front, all fronts with odd k are positive 
fronts and all fronts with even k are negative fronts. Through an analysis similar to the one in the 
previous Section, we may conclude that the dynamics of the fronts can be tracked by the following ODE 
system, similarly to (|3.28[) : 

c(v) if 1 < k < n is odd, (3.38) 
— c(v ) if 1 < k < n is even, (3.39) 
u + (v)L + +u^(l-L + ), L+= ihi+i-hi)- (3-40) 

0<2i<n 

For simplicity of notation, we have dropped the additional subscript showing that the above are leading 
order approximations. As the fronts evolve, it is possible that adjacent fronts will collide. In this case, 
two fronts will disappear, and the dynamics can be continued by renaming the fronts and applying the 
above ODE system with n — 2 fronts instead of n fronts. If front 4>i{t) or cj> n (t) hits either x = or x = 1 
respectively, one can again write down an ODE for the front positions with n — 1 fronts valid after this 
incident. 

As t —> oo in the above ODE system, it is possible that the final configuration will still consist of 
multiple fronts, despite possible annihilations of fronts that may have occurred. At this point, v = v c , 
and all fronts have velocity 0. As far as the intermediate time scale is concerned, these multiple front 
solutions are stable. 

A natural question is whether these multiple front solutions will slowly evolve beyond the intermediate 
time scale. In this long time scale, v is almost exactly equal to v c everywhere. We shall not include a 
detailed analysis of the evolution of multiple fronts in the long time scale, since the analysis and results 
turn out to be very similar to that of the mass-constrained Allen Cahn model, whose long time behavior 
has been studied extensively [HJ HH1 E3- More specifically, the long time behavior of multiple front 
solutions of the equation: 

«e 2 § + /(^ c )-A, A = j f(u,v c )dx (3.41) 
with initial conditions satisfying: 

v c + j udx = K (3.42) 

is the same to leading order to the long time behavior of the multiple front solutions of our system. In 
the mass-constrained Allen Cahn model, multiple front solutions are known to slowly evolve to a single 




d<t>k 



dt 
K = 
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front solution. Thus, multiple front solutions are metastable, and the only genuinely stable solutions are 
the single front solutions. The time scale of this evolution is, however, "exponentially slow" |37j . 

3.4. Higher Dimensions. We have thus far focused our attention on the ID model. We may pose 
similar equations in higher dimensions. Given a bounded domain in Q in K™, we may pose the following 
problem: 

e^ = e 2 Au + f(u,v), (3.43a) 

e-£ = DAv-f(u,v), (3.43b) 

with no-flux boundary conditions on the boundary d£l. If n = 2, we may view this model as modeling 
a top-down view of a cell of (small) uniform thickness. Given that u is a membrane bound species, it 
may be more appropriate in the cell biological context to study the following. Let VL be a smooth domain 
in R 3 whose boundary is the cell membrane T . The membrane bound species u resides entirely on the 
membrane T whereas v diffuses freely inside the cell. We may write down the following equations: 

e— = e 2 A r u + f{u,v) onT, (3.44a) 
e-^ = DAv in 0, (3.44b) 
= f(u,v) on T. (3.44c) 

where Ar is the Laplace-Beltrami operator associated with the surface T and ^ is the normal derivative 
of v taken in the outward direction. It is easily seen that: 

d_ 

dt 



udx + / vdx ) = (3.45) 
Jn 



where J T -dx denotes surface integration over F. 



We now include a brief discussion of the behavior of these systems. We shall first focus on system 
()3.43[) and comment on (|3.44j) later. In the short time scale, the domain Q, will be segregated into a portion 
Q + where u ~ u + and i7_ where u ~ u_. The concentration v is uniform throughout the domain. This 
serves as the initial profile for the intermediate time scale. The dynamics in the intermediate time scale 
may be analyzed analogously to what was done for the ID case. Introduce a stretched coordinate near 
the transition layers, and perform matched asymptotics. This analysis is slightly complicated by the fact 
that the transition layers lie along curves if il C K 2 (or hypersurfaces if il C K n ) and thus, requires the 
introduction of curvilinear coordinates. However, to leading order in the intermediate time scale, effects 
due to the geometry of the transition layers turn out to be higher order. Thus, the fronts advance with 
the normal velocity equal to c(v) where v is spatially uniform and determined by mass conservation. 

The dynamics in the long time scale reduces to the dynamics of the mass-constrained Allen Cahn 
model. The transition layers evolve by mean curvature, with the constraint that the area (or volume) 
of f2+ must be conserved. We refer to [3H1 HI] for further details of the long time behavior of the 
mass-constrained Allen Cahn model. 

The dynamics of p.44[) is similar. In the short time scale, T segregates into domains r + and T_ where 
the concentration of u is u+ and U- respectively, whereas v is uniform throughout S7. In the intermediate 
time scale, the transition curves that reside on T advance so that the normal velocity is equal to c(v). 
In the long time scale, we expect that the transition layers will evolve by geodesic curvature with the 
constraint that the area T + must be conserved. 

Although this long time behavior is interesting from a mathematical point of view, it is not clear 
whether this has any relevance in the cell biological setting since many other biochemical pathways will 
have come into the picture by that time. 

4. Bifurcation Structure of the Wave-Pinning System. As we saw in the previous section, 
wave-pinning behavior occurs for small values e. In this regime, the pinned single front solution (which 
we shall hence forth refer to as the pinned solution or pinned front) is a stable stationary solution of 
the system. A natural question is whether this pinned solution persists as the value of e is increased. If 
e 2 = D (or e = \/~D) in (|2.6[) (the diffusion coefficient of the two species are the same), such stable front 
solutions cannot exist. We thus expect that there is some value of e above which the pinned solution 
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Fig. 4.1. The effect of e on wave shape and existence/stability for reaction terms II2.14H (a-c) and 112, 1611 (d-f) 
with D = 1. (a,d): the shape of the pinned wave with (a) K = 2.8 for e = 0.05,0.1,0.15,0.16 and (d) K = 1.9 for 
e = 0.05,0.1,0.15,0.19. The front gets shallower and broader as e increases, losing stability at (a) e c ft: 0.1621 (d) 
e c ft; 0.1980. We plot the amplitude of u (in b, e) and the mean of v (in c, f) as the pinned solution is continued for 
K = 2.8,2.9,3.0 for (b, c) and K = 1.9,2,2.1 for (e, f). The peaks correspond to saddle-node bifurcations. In (b), the 
amplitude approaches u+(v c ) ft: 1.5150 for the pinned front and decreases as the solution is continued. For K = 2.9, 
the amplitude reaches 0, at which there is a pitchfork bifurcation (see Section In (c), the mean of v is close to 

V = v c = 2.17506 for the pinned solution. In (e), the amplitude is close to 2 for the pinned front. In (f), the mean of v is 
close to 1 for the pinned solution. Note the similarity of (f) with Fig. \4-4\ 



ceases to exist. We simulated (|2.6[) to steady state and gradually increased the value of e. Fig. I4.1f a) and 
(d) show the results of sample computations when (|2.14D and (|2.16[) are used for the reaction term. For 
small e values, there is a gradual change in the wave shape and stall position. Beyond some e c , the pinned 
front disappears and is replaced by a stable spatially homogeneous solution. An interesting feature of 
this transition is that it is "abrupt" : the amplitude of the front (the difference between the maximum 
and minimum values of u) does not vanish gradually as e approaches e c . In Section 031 we shall explore 
the bifurcation structure for (|2.14p and (|2.16[) . In Section |4~21 we focus on obtaining detailed information 
on the bifurcation structure for (|2.16p . In Section 14.31 we indicate other possible bifurcations we may 
expect of the pinned solution. 



4.1. Bifurcation at finite D. For fixed D > and K chosen in a suitable range, there is a stable 
front solution to system (|2.6[) (i.e., the pinned front) for e sufficiently small. As we just saw, there is a 
value e = e c above which this pinned solution can no longer be continued. The value e c is thus a function 
of D and K. 
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To compute e c , we first note that stationary solutions satisfy the following system of equations: 

(4.1a) 
(4.1b) 
(4.1c) 
(4.1d) 

For fixed D and A, we may numerically compute e c by the following procedure. First, set e = eq 
sufficiently small so that the model (|2.6p supports wave-pinning behavior. We discretize (|4.ip and use a 
Newton iteration to find the stationary solution, first for e = eq. The solution to the above equations is 
far from unique, and we thus need a good initial guess to obtain the stationary solution of interest to us, 
the pinned solution. In order to do so, we simulate the time-dependent model (|2.6[) to steady state (i.e., 
run the simulation long enough) with initial conditions suitable to produce a pinned solution as t — > oo. 
This computed steady state serves as our initial guess to find the steady state solution at eq. 

Once we have the pinned solution at eo; we may increase e to e — e\ > eq and use the pinned solution 
at eq as our starting point for the Newton iteration. The stability of the newly obtained stationary 
solution can be monitored by computing the spectrum of the discretization of the linearized operator L 
around the computed stationary solution. This process can be repeated to find e c , the point above which 
a front solution can no longer be obtained. 

Determination of e c using the above procedure, however, is unreliable since the Jacobian matrix for 
the Newton iteration becomes increasingly singular as e approaches e c . We could not obtain values of e c 
beyond an accuracy of approximately 1.0 x 10~ 3 . A well-known remedy for this is to use pseudoarclength 
continuation |32j . Instead of using e as the bifurcation parameter, one uses a pseudoarclength variable 
along the bifurcation curve so that we can successfully continue the solution across a fold singularity. 
This also gives us information on the bifurcation structure at e = e c . 

Results from such a computation are given in Fig. 14.11 (panels b, c for reaction terms (|2.14D and 
panels e, f for reaction terms (|2.16|1 ). For all values of D and A tested, the numerical results indicated a 
fold (saddle- node) bifurcation at e c . The pinned solution is stable until e = e c is reached, and this merges 
with a front solution which has a one-dimensional unstable direction. 

The dependence of e c on D and A is plotted in Figure 14.21 Recall from p.35[) that K m - m = 1 < K < 
3 = A max is the wave-pinning regime for (|2.16j) . For (|2.14[) . the wave-pinning regime is given by (13.3311 
where u± is given in (|2.15[) and v c w 2.17506 is the solution to the equation (see (|2.10[0 : 



u + (v c ) 



1 



I(v c ) = / f(u,v c )du = v c (u + (v c ) - arctan(u + (u c ))) - -(u+(v c )) =0. (4.2) 

Ju-(v c ) 

The wave-pinning regime is K m i n < K < K max where K m i n = v c ,-Kma X ~ 3.6901. For both (|2.14[) and 
pT6]l . we sampled K between (3K nlin + if max )/4 <K < (K min + 3if max )/4. 

For both reaction terms, if we fix D, we see that there is a value K = K m at which e c reaches a sharp 
peak. The value of K m is located roughly at (K min + A" max )/2. This may be rationalized as follows. 
The asymptotic calculations were based on the front being sharp and away from the boundaries. When 
K = (A' lmn + A max )/2, the front is positioned approximately in the middle of the domain, and thus 
farthest from the boundaries. One may thus argue that A = (A' m i n + A max )/2 would tend to "maximize" 
the range of validity of the asymptotic calculations. The reason for the peaked appearance of the e c plot 
for fixed D will be addressed toward the end of the next section. 

The computed results indicate that e c is uniformly bounded in D and A. In particular, we observe 
that, for fixed A, e c (D,K) tends to some value as D is taken large. This serves as one motivation for 
studying the limit D — > oo. 

4.2. Bifurcation Diagram in the limit D — > oo. Here, we study the bifurcation structure of the 
following system: 

du ^d 2 u ,, . . 

(4 ' 3a) 

v = K- udx, (4.3b) 
Jo 
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Fig. 4.2. Composite fam- 
ilies of two-parameter bifurca- 
tion diagrams showing the wave- 
pinning regimes (always below 
the displayed curve(s)) in the 
Ke and De planes for kinetics 



lITTB (in a, b) and STTEl (in 

c, d). In (a) and (c), the criti- 
cal value of t, e c , is plotted as 
a function of K for fixed D. 
The five solid lines correspond 
to D = 0.25,0.5,1,2,4. The 
dashed line in (c) is the D — > 00 
curve (computed separately, see 
Section \4-8ty ■ In (b) and (d), e c 
is plotted as a function of D for 
fixed K where K = 2.8, 2.9, 3.0 
in (b) and K = 1.9,2.0,2.1 
in (d). The D-axes are scaled 
logarithmically. 



0.15 




0.25 



0.15 
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2 2.2 2.4 
K 




where u satisfies no-flux boundary conditions at x = and x = 1. The above system should be seen as 
the limiting system of (|2.6[) when D — > 00. In this limit, v is spatially uniform and hence the evolution of 
v is determined completely by the mass constraint. By retracing the arguments of Section [3.21 the reader 
can easily check that the above system exhibits the same wave-pinning behavior as (|2.6[) . In fact, the 
analysis is somewhat simpler. When D is finite, we deduce from our asymptotic ansatz that v is spatially 
uniform to leading order. In the case of (|4.3[) . the spatial uniformity of v is automatic. We note that 
system (|4.3p is often referred to as the shadow system of (|2.6p , and that it has been shown to provide 
insight into the behavior of the original system [551 SI IZ] • 

In studying the bifurcation of system (|4.3[) . we may proceed similarly to the finite D case of the 
previous section. However, we shall take a different approach that will allow us to obtain a far more 
detailed picture of the bifurcation structure. 

The steady state solutions to (|4.3[) satisfy: 



= , 2 g + /(«,„), 
du 

— = at x = 0, 1, 
ox 



v = K 



udx. 



(4.4a) 
(4.4b) 
(4.4c) 



We shall view the first equation as an ODE for u to be solved in the "time" variable x. It is slightly more 
convenient to use r = x/e as our "time" variable. Rewrite the first equation as a system of first order 
ODEs: 



u T = w, 

W T = ~f(u,v). 



(4.5a) 
(4.5b) 



Note that this ODE system possesses an "energy". Multiplying both sides of (|6.3[) by w = u T and 
integrating, we obtain 



= F(u,v,B), F(u,v,B) = -B + F Q (u,v), F (u,v) = 2f(s,v)ds 



(4.6) 



where B is an integration constant. Consider the u — w phase plane that corresponds to system (|4.5[) . The 
solution curves of (|4.5[) coincide with the level curves of the function w 2 = F (u, v) (see Fig. I4.3[) . Recall 
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Fig. 4.3. Typical shapes 
of the functions y = Fq(u,v) 
(panels a,c) and level curves of 
w 2 = Fo(u,v) (panels b, d) in 
the uw phase planes for kinetics 
H2.16II for v = 1 ( right ) and v = 
9/8 (left). It is clear that there 
can only be a half-loop trajec- 
tory when Fq(u, v) = B has four 
distinct solutions. This happens 
when B min < B < B max . Note 
that, as B ~\ B m i n , the half 
loop approaches either a hetero- 
clinic or (half of) a homoclinic 
orbit. As B /" B max , the half 
loop approaches the neutral cen- 
ter (u, w) = (u m ,0). 



that the function f(u, v) is bistable in u for fixed v satisfying v m i n < v < i> max (the bistability condition). 
We shall be interested only in analyzing cases in which v falls within this bistable range. In this case, 
the function y = Fq(u, v) for fixed v has the form of a double well potential, whose local minima are at 
u = and whose local maximum is at u = u m . A stationary solution of system (|2.6p satisfying 

equation (|4.4b|) corresponds to a solution trajectory in the u — w phase plane that starts and ends at 
the it-axis (or w = u T = eu x = 0). It is clear that there can only be such a trajectory if B = Fq(u,v) 
as an equation for u has four distinct solutions (see Fig. I4.3[) . Let the two middle roots be uq and u\ 
(we assume uq < u\). Then, stationary single front solutions correspond to the "half loop" trajectory 
that connects (it, w) = (uq, 0) and (u±, 0). Wc sec immediately that such stationary single front solutions 
must be either monotone increasing or decreasing. In fact, the only stationary solutions (|2.6[) can have 
are multiple "half loop" trajectories that correspond to periodic multiple front solutions. 

For every (v, B) such that F(u, v, B) = has four solutions in u, we have a corresponding half loop 
trajectory. We have, therefore, a two parameter family of half loop trajectories, and hence of possible 
stationary single front solutions. In the following, wc shall simply refer to stationary single front solutions 
as front solutions. Let T> v b be the range of (v,B) values for which F(u,v,B) = has four solutions. It 
is clear that (see Fig. 14. 3|) : 

T^vB = {(v, B) G K 2 |u m in < V < v max , B min (v) < B < B mllx (v)} ^ ^ 

B min (v) = maxF (u±(t;),w), B max (v) = F (u m {v),v) 

The expression for B max denotes the greater of the values Fq(u + (v), v) and Fq(u- (v), v). We thus have a 
correspondence between (v, B) values in T> v b and half-loop trajectories. The set V v b exhausts all possible 
half-loop trajectories but one cannot say in general whether this correspondence is one-to-one. However, 
we do expect this to be true if / is not too pathological. For reaction terms (|2.16l) . (|2.17[) or (12.71) . it is 
not difficult to see that this is indeed the case. We shall henceforth assume that this correspondence is 
one-to-one. 

The task of finding front solutions has now been reduced to finding the suitable half-loop trajectories 
that satisfy (|4.4b|) and (|4.4cp . First, consider (|4.4b|) . Half- loop solutions automatically satisfy w = u T = 
and hence u x = at the endpoints, but this does not necessarily mean that the endpoints correspond 
to x = 0(t = 0) and x = l(r = 1/e). We must thus impose the condition that the domain length is 
1. Suppose the front solution has value uq at x = and u\ at x = 1, uq < u\ (we shall henceforth 
assume that our front solution is always monotone increasing, unless noted otherwise). The domain 
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length condition reduces to: 

1 f Ul dx f Ul du 

dx= ^du = e J U ee el(v, B) (4.8) 

Juo du J U0 y/F(u,v,B) 

where we used (|4.6[) . The above change of variables is valid because we know that the stationary front 
solution is monotone increasing. Note that uo and u\, being the middle roots of the equation F(u, v, B) = 
0, are functions of v and B. Hence, the above integral is a function of v and B. Condition (|4.4c[) can, 
likewise, be written in the following form. 

f Ul udu 

K = v + f m m = v + eJ(u ' B) - (4 - 9) 

Ju„ ^F(u,v,B) 

We have thus reduced (|4.4j) to the two integral constraints (|4.8[) and (|4.9p . Furthermore, the integral 
constraints incorporate the fact that we seek single-front solutions; (|4.4p is satisfied by any stationary 
solution. Given e and the total mass K, we may solve (|4.8[) and (|4.9D for v and B, which in turn uniquely 
determine the half-loop trajectory, and hence, the solution u. We note that a similar reduction is possible 
even when D is finite. We describe this in Appendix 16.11 When D is finite, however, it seems somewhat 
difficult to use this reduction to great effect. 

Conditions (|4.8[) and (|4.9p may be used in place of (|4.4[) as the basis for continuing the pinned 
solution. The use of conditions (|4.8p and (|4.9p has a computational advantage over the direct use of 
(|4.4p since the former is a much smaller system to solve than the latter. We do note, however, that 
the numerical evaluation of the integrals I(v 1 B) and J(v, B) is not entirely trivial, especially when B is 
close to i? m in- This is related to the fact that the half- loop trajectories come very close to heteroclinic or 
homoclinic orbits on the u — w phase plane. The techniques used to overcome this difficulty are discussed 
in [TT]- 

A more interesting use of the above conditions is the following. Since e ^ 0, we may eliminate e from 
(GLl and (HU. We have: 



>K 



(v, B) = (K - v)I{v, B) - J{v, B) = 0. (4.10) 



If we can find the zero set of Qk{v,B) where {v,B) € T> v b, we will have obtained all single front 
stationary solutions for a fixed value of K (with v in the bistable range) regardless of whether it arises 
as a continuation of the wave-pinned solution. 

Any point on this zero set corresponds to a different front solution, and the value of e can be recovered 
by using (|4.8p . Consider the map: 

M : (v,B) — > (M,e) = (M(v, B), (I(v, B))^ 1 ) (4.11) 

where the function M(v, B) is chosen so that the map Ai defines a homcomorphism on T> v b- Note that 
the choice of M is far from unique; we shall see that M(v,B) = v works well for (|2 . 16|) and (|2.17l) . 
Half-loop trajectories can then be parametrized by (M, e) instead of (w, B). The zero-set of Qk(v, B) in 
T> v b can be mapped by Ai in a one-to-one fashion to yield a bifurcation curve on the M — e plane. 

Up to now, the treatment has been fully general. We now apply this methodology to the case when 
the reaction term is given by (|2.16p . We shall be interested in obtaining the bifurcation diagram when 
1 < K < 3, the wave-pinning regime (see (|3.35p ). First, we note that < v is the bistable range. The 
domain T> v b is therefore an unbounded set, making it difficult to uniformly sample points in T> v b to 
determine the zero set of Qk(v, B) and hence the bifurcation curves. The following considerations allow 
us to restrict our search to a much smaller set. Given that uq and u\ are the two middle roots of equation 
F(u, v, B) = 0, we have: 

= U-{v) <u < u rn (v) = 1 < m < u + (v) = 1 + v. (4-12) 

Therefore, 

< u < I udx < ui < 1 + v. (4-13) 

Using (|4.1dp . 



o 



v < v + / udx = K < 1 + 2v. (4-14) 
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Therefore, we may restrict our search of the zero set of Qk(v,B) to the following range: 

<v< K. (4.15) 

We note in passing that a similar argument can be used to show that any single front stationary solution 
to (|4.3[) without restriction on v (for 1 < K < 3) must in fact satisfy v > and hence (|4.15[) . 

We thus numerically evaluate Qk(v, B) at sample points in the range T>^ B = V v b H {^-^- < v < K] 
to find the zero set of Qk{v, B). More specifically, we fix v and sample B uniformly within the admissible 
range. If there are adjacent sample B points for which Qk(v, B) changes sign, a zero is obtained between 
these values by bisection. This procedure is repeated for v values uniformly sampled in f|4. 15[) . Where 
the zero set has a complicated structure, sampling is refined to clarify this structure. Once the zero-set 
is obtained, we use the map Ai with M(v,B) = v (see (|4. 1 1 [) ) to obtain a bifurcation curve in the v — e 
plane. Computational evidence indicates that e = (I(v, -B)) -1 is an increasing function of B for fixed v, 
and thus M. : (v, B) i — > (v, e) is a homeomorphism on D v b- 

We can explicitly obtain the region T> ve = M(D v b) by studying the integral I(v, B). Assuming that 
I(v, B) > is a decreasing function of B for fixed v, we have only to know the limiting values of I(v, B) 
as B — > B m i n (v) and B meix (v) (see (|4.7p ). As B — > B min for fixed v, the half loop trajectories approach 
(half of) a homoclinic orbit or a heteroclinic orbit in the u — w plane (see Fig. I4.3|) . In either case, the 
total "time" it takes for the orbit to complete the half loop increases as B — > B m - m . Thus, / — > oo as 
B — > B m \ n . On the other hand, when B — > i? max , the half-loop trajectory approaches the neutral center 
(it, w) = (u m , 0) = (1, 0) in the u — w phase plane. We may easily compute: 

lint I(v,B) = ^=. (4.16) 

Therefore: 

V m = {(v, e) G M 2 |0 < v, < e < V^/n}. (4.17) 

As one approaches the parabolic edge of T> ve , the amplitude(= ui — uq) of the front solution tends to 
and approaches the spatially homogeneous steady state (u, v) = (l,u). In fact, the parabolic edge is the 
only place where the amplitude tends to in T> V€ . Combining this with (|4.17p with (|4.15p . we obtain an 
upper bound e < yKj-K for the existence of single front solutions. As we shall see, this bound is not 
sharp. 

To study the stability of the stationary solutions corresponding to points on the zero set, we must 
compute u explicitly. Once we know v and B, we can find uq and u\. We can then numerically solve the 
initial value problem (|4.5[) with the initial values u(0) = uq and w(0) = 0. Up to numerical error, the 
computed solution should, by design, satisfy u{t = 1/e) = u(x = 1) = u\ and w(l/e) = u t (t = 1/e) = 
eu x (x = 1) = 0. We can then linearize about u the operator on the right hand side of (|4.3[) . By examining 
the spectrum of (the discretization of) this linearized operator, we can determine linear stability of the 
steady state u. 

The resulting bifurcation curves on the v — e plane are given in Fig. 14.41 When Ky^2,we found that 
there was at most one front solution that corresponds to each value of v (or equivalently, Qk(v,B) = 
had at most one solution in B for fixed v). We first consider the case K < 2 (left panel). For 
small values of e, there are three front solutions. In order of increasing v, we denote these solutions 
(it, v) = (u~(x),v~), (u™ p (x),vj p ) and (u^(x),v^), which we call the minus, middle and plus branches 
respectively. 

The pinned solution corresponds to the middle branch (it™ p (x), v ™ p ). The value v^ p approaches 1 
as e — > 0. We know from our asymptotic calculations that the integral J 1 f(u,v)du vanishes to leading 
order when the wave stalls. This happens when the three roots u±(v) and u m (v) are equally spaced, 
which corresponds to v = 1. In the u — w phase plane, u™ p approaches a heteroclinic orbit that connects 
the two saddle points (u,w) = (0,0) and (u, w) = (2,0). 

The other two front solutions are unstable and have a one-dimensional unstable direction. The values 
v~ and w+ approach v = (K — l)/2 and v = K respectively as e — > 0. Let us consider the plus branch. 
As e — > 0, B approaches B m i n . In the u — w phase plane, uf approaches (half of) a homoclinic orbit 
that originates and ends at the saddle point (u, w) = (0,0). As uf approaches this homoclinic orbit, the 
amount of "time" that the solution stays close to the saddle point increases, so that uf is very close to 
for much of the interval < x < 1. Near x = 1, there is a sharp transition zone in which uj makes 
a steep increase to u±. This transition zone becomes increasingly narrow as e — > 0. We may say that 
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Fig. 4.4. Bifurcation diagrams for cubic kinetics II2.16I in the v — e plane for K < 2, K = 2 and K > 2 (from left 
to right). When K ^ 2, iae middle branch is stable and the others are unstable, except for a small region 2 < K < fa 
2.00672 (details in text and Fig. |^.5| ). W/ien X = 2, middle and minus branches meet at a transcritical bifurcation (TC, 
[vtc, e*c) ~ (1, 0.23530) J and exchange stability. The values and tend to (K — l)/2 and if, respectively, as e — > 0. 
27ie slopes of the line of tangency of the curves (v^~,e) and (v^~,e) as t — > are calculated in 114,1911 and 14.1811 . The 
pitchfork bifurcation (PF) occurs at (v p f,e p f) = (K — 1, sJK. — l/7r). JVcrfe iaat tfte above bifurcation diagram only shows 
the monotone increasing front solution. At the pitchfork bifurcation, this meets with the monotone decreasing front as well 
as the spatially homogeneous state. SN: saddle-node bifurcation. 

the solution (uf(x),v+) approaches the stable homogeneous steady state (u,v) = {0,K) as e — > 0. The 
convergence of uf(x) to is only uniform outside of an arbitrarily small neighborhood of x = 1. We can 
use the above phase plane information together with (|4.9I) to find the following approximate expression 
for when e is small: 

w c + = K - a + e + o(e), 

a + = |m J(K,B) = 2^1n (^=^) , (4.18) 
At = I (2(2 + K) ± ^2(2K + 1)(K -1 

where o(e) is the usual Landau symbol denoting an expression that tends to as e — > 0. Note that 
the expression inside the square root in f3± is positive since K > 1. The validity of this expression is 
supported by computational results. 

The situation for the minus branch is similar. As e — > 0, v~ — > (K — l)/2. In the u — w phase plane, 
u~ approaches half of the homoclinic orbit originating from the saddle point (w, w) = {{K + 1) /2, 0). The 
value of u~ is close to (A" + l)/2 for most of < x < 1 except for a small neighborhood around x = 0. 
The function u~ converges uniformly to [K + l)/2 on any set outside an arbitrarily small neighborhood 
around x = 0. Similarly to (|4.18[) . we can obtain the following expression for v~: 



A - 1 
2 



a_e + o(e). 



^ i °(^^)'-=5r ± v5« 3+A -» (3 - A - ) 



(4.19) 



Note that the expression inside the square root in j± is positive since A < 3. The validity of the above 
expression is supported by numerical results. 

As e is increased, there is a value e = at which the middle and plus branches meet in a saddle-node 
bifurcation. At this point, the disappearance of the stable front solution (the middle branch) is "abrupt" 
in the sense that the amplitude of u is non-zero as the bifurcation is approached. This can be seen from 
the fact that this saddle- node bifurcation occurs in the interior of T> vc (see discussion after (|4.17[) ). 

The minus branch can be continued until it merges with a spatially homogeneous unstable steady 
state. This happens at the parabolic edge of T) ve (see (|4.17p ). Here, there is a pitchfork bifurcation 
at which the homogeneous steady state gives rise to two unstable front solutions, one that is monotone 
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increasing and the other monotone decreasing. Note that, in Fig. 14.41 only the bifurcation diagram 
of the monotone increasing front solution is plotted. There is an identical bifurcation diagram for the 
monotone decreasing front, and these two solutions meet with a spatially homogeneous steady state at a 
pitchfork bifurcation. This homogeneous steady state corresponds to (it, v) = (1,K — 1). Using (|4.17D . 
the corresponding e value e p f is equal to: 



Linearize the operator on the right hand side of (|4.3[) around this homogeneous steady state, and call 
this operator C We can also obtain the above value by considering the spectrum of C. At e = e p f, one 
of the eigenvalues of C corresponding to the wave number k = tt becomes positive (see Section 13.11 in 
particular, equation p.lOj) in the D — > oo limit). For 1 < K < 3, expression (|4.20p gave the least upper 
bound of the range of e for which (not necessarily stable) single front solutions exist. 
We now turn to the case K = 2. Let us first take a look at (|4.10[) . If v = 1, we have: 



The function F(u, 1,B) is symmetric about u = 1 and thus the same is true for uq and u\ (i.e. (uq + 
ui)/ 2 = !)■ The above integral is therefore equal to whenever it is well-defined. Therefore, all points 
such that v = 1 in T> vc (i.e., (v, e) = (1, e), < e < 1/tt) are part of the bifurcation curve for K = 2. For 
e small, this v = 1 branch of solutions corresponds to the pinned front, which is thus stable for small e. 
We shall denote this branch by (u™ p ,vj p ) and call this the middle branch. As e / 1/tt, we expect the 
middle branch to merge with an unstable homogeneous steady state at a pitchfork bifurcation, just like 
u~ for K < 2. An unstable solution cannot give rise to two stable solutions in a pitchfork bifurcation. 
We must conclude that the middle branch is unstable when e is close to 1/tt. This suggests that there 
must be an intermediate e value between and 1/tt at which the middle branch loses stability. This is 
indeed the case. 

Just as in the K < 2 case, there are three single front solutions in the K = 2 case for small e. We 
shall refer to them in the same way as in the K < 2 case. As we saw, v™ p is always equal to 1. As e is 
increased, the minus and middle branches meet in a transcritical bifurcation at e = e tc ~ 0.2353. Above 
etc, the minus branch becomes stable and the middle branch loses stability. At e = e+ w 0.2419, the 
minus and plus branches meet in a saddle-node bifurcation. 

The K > 2 case is similar to the K < 2 case except for some fine details. When e is small, we have 
three front solutions which we name in the same fashion as in the K < 2 cases. The middle branch 
merges with the minus branch at e = in a saddle-node bifurcation. The branch plus branch merges 
with the spatially homogeneous solution at e p f = \JK — 1/tt (sec (|4.20[0 in a pitchfork bifurcation. 

An interesting detail in the K > 2 case is that there is a small window 2 < K < K p w 2.00672 for 
which the plus branch has a stable portion (see Fig. 14. 5[) . The existence of such a portion is implied 
by the structure of the bifurcation diagram at K = 2. The saddle-node bifurcation at e+, should persist 
beyond K = 2 since saddle-node bifurcations are robust under perturbations. On the other hand, when a 
transcritical bifurcation is perturbed, it will generally give rise to zero or two saddle node bifurcations (see 
[TT] for example). When K = 2 is perturbed to K < 2, the transcritical bifurcation does not give rise to 
any saddle node bifurcations. If perturbed to K > 2, it gives rise to two saddle node bifurcations, one of 
which occurs at e = e^. We shall name the other e value e = e® n . Both bifurcation points corresponding 
to e = 6g n and e^ n lie on the (uf,v^) branch. 

For 2 < K < K p , there are three single front solutions over the range e® n < e < e^ n , which we refer 
to as the (+,l),(+,2) and (+,3) branches respectively in order of increasing v. The (+,1) and (+,2) 
branches meet in a saddle-node bifurcation at e® n and (+, 2) and (+, 3) branches meet at e+,. The (+, 1) 
and (+, 3) branches are unstable whereas (+, 2) branch is stable. For 2 < K < K p , therefore, there is a 
small window of e values for which there is a stable front solution that cannot be reached by continuing 
the pinned front solution. For K > K p , the plus branch does not have a stable portion. At K = K p , the 
saddle-node bifurcation points merge and disappear. 

In Fig. 14.61 we show the (K, e) parameter region in which there is a stable single front solution. This 
should be seen as a refinement of the e c plot in Fig. 14.21 that we obtained for finite D. The region is 
peaked at approximately K = 2, but with some fine structure coming from the small window of front 
solutions that exist for 2 < K < K p . The peaked geometry of this region comes from the fact that the 
saddle-node bifurcations at which the pinned solution loses stability are different for K > 2 and K < 2. 
At K = 2. we have a transcritical bifurcation that separates these two regimes. 




(4.21) 
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Fig. 4.5. Bifurcation dia- 
gram for cubic kinetics 112. 161 1 
for 2 < K < K p ft 2.00672 
(K = 2.001 in this example) 
plotted in the v — e plane. The 
full diagram is on the left panel, 
a part of which is magnified on 
the right panel. The (+,1) and 
(+, 2) branches (labeled uj"' 
and v^~' 2 respectively) come to- 
gether at the saddle-node bi- 
furcation denoted SNO and the 
branches for (+, 2) and (+, 3) 
branches (labeled v^' and v^' 
respectively) come together at 
SN+. The (+,2) branch is sta- 
ble. At SN-, middle and minus 
branches come together. The 
dotted lines are the bifurcation 
curves at K = 2. 



Fig. 4.6. Two parameter 
bifurcation plots for cubic kinet- 
ics 112.161 1. A stable front solu- 
tion exists for parameter values 
in the K — e parameter region 
bounded by the curve and the K- 
axis (left panel). On the right 
panel, the "tip" of the curve 
is magnified. At (K p ,e p ) ft 
(2.00672, 0.24474) the values ef„ 
and e® n come together. At 
(2, et c ), etc ft 0.23250, e~ n and 
e® n come together at the trans- 
critical bifurcation point. 

It is not clear how much of the insights we obtained for D — > oo can be carried over to the finite D 
case or to the reaction term (|2.14[) . It seems plausible, however, that much of what we learned does indeed 
carry over. For example, we expect that there is a K value (that depends on D) at which the pinned 
solution undergoes a transcritical bifurcation (rather than a saddle-node bifurcation) as e is increased. 
The peaked appearance of the e c plot of Fig. 14.21 serves as circumstantial evidence for this claim. 

4.3. Other Possible Bifurcation Structures. We now have a clear picture of the bifurcation 
structure for reaction term (|2.16[) . especially when D — >• oo. Given the broad similarity of the e c plots for 
(|2.14[) and (|2.16[) (sec Figure |4~2")) . it is natural to expect (|2.14[) to also have a bifurcation structure with 
features similar to (|2.16[) . This raises the question of how general our findings are. For other reaction 
terms that support wave-pinning, there is no reason to expect the full bifurcation structure to be similar. 
In particular, we can raise the following question. Except at K = 2, the pinned front was seen to undergo 
a saddle-node bifurcation in the case of (|2.16p . D — > oo. This bifurcation was "abrupt" in the sense that 
the front amplitude tends to a non-zero value as the bifurcation point is approached. Is the saddle-node 
bifurcation the only generic way in which the pinned front is lost? In particular, is it generically the 
case that the bifurcation is "abrupt"? The answer to both questions turn out to be negative. We shall 
demonstrate this with a description of the bifurcation structure for the reaction term (|2.17[) . We shall 
see that the pinned front can arise generically via a pitchfork bifurcation from a spatially homogeneous 
state. The exposition will be kept brief since much of the analysis proceeds along lines similar to that 
of the previous section. We shall only discuss the D — > oo case. The case of D finite is expected to be 
similar. We note in particular that computational examples can be produced in which such bifurcations 
occur for finite D. 

As we saw in Section 13.11 an interesting feature of the reaction term (|2.17[) is that the middle 
homogeneous steady state (u m ,v) can be stable. A similar conclusion is true in the D — > oo case. This 
happens when a > 1 in (|2.17p . We shall concentrate on this case. When a < 1, the full bifurcation 
diagram turns out to be quite similar to that of (|2.16[) (the generic bifurcation is the saddle-node). We 
shall not discuss this case here. 

As can be easily checked, — oo < v < oo is the bistable range, and — 1 < K < 1 is the range for which 
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Fig. 4.7. Bifurcation diagrams for the reaction term i2.1H in the v — e with a = 2 and indicated values of K. Left: 
K = 0.1 < K r ss 0.19498, Middle: K r < K = 0.22 < K q ss 0.22510, v v J is represented by the small portion of the curve 
between the pitchfork bifurcation (PF) and the saddle-node bifurcation (SN). Right: K = 0.3 > K q . 



wave-pinning can occur. We focus on these values of K. 

Let us first study the spatially homogeneous steady states of the system for fixed K. Given v, u 
must be either u = u +1 u- or u m . There is one spatially homogeneous steady state each for u_ and u+: 
(u_, v) = (—1, K + 1) and (u + , v) = (1,K — 1). Let us consider the spatially homogeneous steady states 
that correspond to u = u m . Given that the total mass must equal K, we have the following equality: 

v + u m (v)=v- aV =K. (4.22) 
y/1 + [avY 

It turns out that this equation can have three solutions in v for fixed K if a > 1 and K satisfies: 

-K q <K < K q , K q = i(a 2 / 3 - If' 2 . (4.23) 

It is clear that K q is always smaller than 1. Let us call these three solutions i>~ < w° n < u+. We may 
adapt the calculations of Section 13.11 to the D — > oo case. It can be checked that to = f u — f v < at 
(it^, i>^J = (u m (u^J, v^), and therefore, that this is a stable steady state so long as: 



e > = e° f (4.24) 

Note that the above expression can be obtained by taking D — > oo in p.lOp . We name the right hand 
side £p f in anticipation of our results to be discussed below. The other two homogeneous states are always 
unstable. When \K\ > K q (|4.22[) has only one solution. 

The bifurcation diagram in the D — ¥ oo limit can be obtained in a procedure similar to the treatment 
of (|2.16l) in the previous section. The possible bifurcation diagrams in the v — e plane are given in Fig. 
14.71 where we have taken a = 2 in (|2.17[) . Only the case K > is shown. Given the symmetry of the 
reaction term (|2.17[) . the bifurcation diagram for —K can be obtained by flipping the bifurcation diagram 
for K about the e axis. 

For all values of — 1 < K < 1, there are three single front solutions when e is sufficiently small. Just 
as in the previous section, we shall denote them by (u~ , v~), (uj p , vj p ), (u+, u+) and refer to them as the 
minus, middle and plus branches. There is a constant < K r < K q (that depends on a) such that, when 
< K < K r the middle branch merges with the stable homogeneous solution (u® n , w° ra ) in a pitchfork 
bifurcation. This happens at e = ep f whose analytical expression was given in (|4.24[) . The pinned front 
solution is stable up to this pitchfork bifurcation. Note that this is only possible since (w^j, uJJJ is a stable 
steady state for e > ep f . The minus and plus branches are unstable and merge in pitchfork bifurcations, 
respectively, with the unstable homogeneous states (u^j,^). 

When K r < K < K q , the situation for the plus and minus branches does not change. However, the 
middle branch now loses stability in a saddle-node bifurcation with the solution branch (u pt , v pi ) that 
arises from a pitchfork bifurcation from the homogeneous state (m^,^). This branch is unstable. The 
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difference between < K < K r and K r < K < K q is whether the pitchfork bifurcation at is 
subcritical or supercritical (see Fig. I4.7|) . In fact, we encountered a similar bifurcation for (|2.14[) when 
D = 1 and K = 2.9 (sec Fig. O(b) and (c)). 

For K > K q , the middle branch loses stability in a saddle- node bifurcation with the minus branch. 
The plus branch merges with the unstable homogeneous solution (u+,v+). The case K = K q is highly 
degenerate and atypical, and we thus omit the details here. 

Assuming that the above bifurcation picture is valid for all values of a > 1 (an observation supported 
by computational evidence), we can compute K r as the value of K at which the pitchfork bifurcation 
at (u^j,^) changes from being subcritical to supercritical. We can then obtain an explicit analytical 
expression for K r , whose derivation we defer to Appendix 16.21 We note that values of K r obtained by 
bifurcation computations match perfectly with the analytical expression we now state. Consider the 
equation: 



5 a 3 
3 + s 3 - a + 8(s 2 - 1) 



0. (4.25) 



It can be shown that there is just one solution to the above equation in the range 1 < s < ^fa. Take this 
root and let: 



J s 2 — 1 av r 

-, K r = v r -. (4.26) 



a s 



The value v r is the value of when K = K r . We thus have an expression for K r as a function of a. Wc 
can see that K r — > 1 as a — > oo. Letting a = sf-^a, we can rewrite (|4.25p as: 



31a 3 - 7 = a2/3 _ 



cr 2 (40ct 3 - 16) 



(4.27) 



We see that a ^16/40 = -^2/5 as a ^ oo. Using (jT^6|) . 



y? - ! \/s 2_ ~T \ / Va 2 / 3 CT 2 - 1 /a 2 / 3 cr 2 -l 

hm + \ — — 

a->oo \ a V a z ' a a 



K r = I — + — I = ( : + \l o/ 3fT 2 I = L ( 4 -28) 



In other words, the range of K over which the pinned solution merges with a stable homogeneous solution 
increases with a, covering the entire wave-pinning regime (—1 < K < 1) as a — > oo. 

For reaction term (|2.16[) . the only generic bifurcation through which the pinned solution is lost was of 
saddle-node type. Reaction term (|2.17|) is an example in which the pinned solution can be generically lost 
by merging with a stable spatially homogeneous state. As a — > oo, this is the case for most values of K in 
the wave-pinning regime. Although these examples give us interesting insight into the possible bifurcation 
structure of (|2.6[) . it is difficult to draw conclusions that may be applicable to arbitrary reaction terms. 
Our study in the present section suggests a general connection between the stability of homogeneous 
states of type (u m (v), v) and the type of bifurcation at which the pinned solution is lost. 



5. Discussion. Previously we have studied the reaction-diffusion model (|2.ip with kinetics (|2.3[) . 
motivated by an investigation of the redistribution of polarity proteins (Rho family GTPases) in eukaryotic 
cells. These switch-like proteins interconvert between an active and an inactive form and diffuse across 
the cell. The appearance of a small parameter, e in this problem stems from the membrane confinement 
of one of the species (the active form) , which tends to reduce its rate of diffusion by orders of magnitude 
relative to the other form. The inactive form diffuses rapidly, i.e. D = 0(1). Conservation of total amount 
of protein (K to tai, and in dimensionless form K) stems from the fact that there is no net production nor 
loss of total protein on the timcscalc of interest. 

In previous work on this biological problem, we had postulated bistability based on positive feedback 
between the presence of the active form and its own activation. This led us to find a phenomenon of wave- 
pinning, which could account for polarization in response to large enough stimuli }21) . We found that the 
phenomenon depends on the ratio between the two diffusion coefficients being small enough. Many other 
proposed models for cell polarization are based on diffusion-driven, Turing- type instabilities [3H1 [Ml US] > 
in which a state that is stable in the absence of diffusion is destabilized in its presence, a mechanism 
fundamentally different from the wave-pinning mechanism considered here. Our main motivation has 
been to understand this phenomenon from a mathematical point of view. 

We first analyzed wave-pinning exploiting the smallness of e using matched asymptotic analysis. We 
identified three key properties the reaction term must satisfy (bistability, homogeneous stability and the 
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velocity sign conditions, see Section [2]) in order for the system to exhibit wave-pinning. Both (|2.3[) . as 
well as the simpler (|2.16|) satisfies these properties and thus supports wave-pinning. The analysis allowed 
us to determine the range of K values for which wave-pinning is possible. Furthermore, we were able 
to reduce the RD system to a simple differential algebraic system for the front position, whose explicit 
form could be calculated in the case of (|2.16D thanks to its algebraic simplicity. This reduction gives 
an excellent approximation of the original system as e is made small (Fig. 13. ip . We briefly discussed 
the long-time behavior of our system as well as its higher dimensional generalizations. We argued that 
the long-time behavior is analogous to that of the mass-constrained Allen-Cahn model, whose properties 
have been well-characterized [HI S3 EM EH] • 

As e is increased, the matched asymptotic calculations are no longer valid, and the pinned front is 
eventually lost. This led us to examine the bifurcation structure of the system. For finite -D, we did 
so using pseudoarclength continuation on the full PDE system (Figs. 14. 1114. 2p . Reaction terms (|2.14j) 
and (|2.16[) revealed a similar bifurcation structure. We found that the pinned front was always lost in a 
saddle-node (fold) bifurcation, and delineated the parameter region in the K — e plane for which wave- 
pinning was possible (Fig. 14. 2[) . We obtained a complete bifurcation picture for single front solutions 
in the limit D — > oo for the reaction term (|2.16[) . using a method related to the "time map" technique 
[34ll6]. We found that there is a transcritical bifurcation for a particular value of K{= 2) (Fig. 14. 4p . This 
value acts as a "watershed" explaining the cusp-like form seen in Fig. 14.21 Other bifurcation pictures 
are possible. In the case of (|2.17p . as shown in Fig. 14.71 the pinned front solution can be lost through 
a pitchfork bifurcation. The possibility of such a bifurcation depends on the stability of the "middle" 
homogeneous steady states. It seems to be difficult to give a general account of the bifurcation structure 
for wave-pinning systems. We hope the two scenarios we identified are representative of what can be 
expected. 

The simplicity of our model and the universality of reaction-diffusion systems in biology, chemistry, 
and physical settings suggests that such wave-pinning phenomena may be quite ubiquitous |19U311l35ll42j . 
In this paper, our motivation stems from cell polarization and the biochemistry of Rho proteins, and the 
variables u and v correspond to active and inactive forms of one Rho protein. More detailed models for 
the dynamics of these proteins, with mutual interactions and effects on the actin cytoskeleton [T2] [T51 [5] 
show similar wave-pinning phenomena, but their complexity makes a full mathematical analysis much 
harder. 

We conclude with a discussion of possible biological implications. The small parameter e exploited 
in our analysis depends on several biological parameters including rates of diffusion D u , reaction rj, and 
domain size L. The necessary condition e <C 1,D ps O(l) is satisfied by virtue of the large difference in 
diffusion of the membrane-bound active Rho protein and inactive form that diffuses freely in the cytosol. 
Normally, these rates of diffusion differ by 100-fold. Assuming a typical cell diameter of 10 /im, reaction 
timescale ?/ = Is -1 , and diffusion coefficients D u = 0.1/im 2 s -1 and D v = 10/im 2 s -1 , the dimensionless 
constants are e ~ 0.03 and D sa 0.1. This is within the wave-pinning regime for the Hill function kinetics 
(|2.3p . However, increasing the diffusion coefficient of the active form tenfold to D u = 1 /im 2 s _1 , or slowing 
down the rate of interconversion 77 to 0.1 s _1 , or decreasing the cell size to L ps 3^m leads to e w 0.1 
and Dsil, putting the Hill function kinetics system into the bifurcation regime where wave-pinning and 
hence polarization would be lost. 

Such predictions are experimentally testable. Cell fragments capable of polarization [35] could be 
made successively smaller to test the effect of domain size. Manipulating the amount of Rho protein 
could test the predicted effect on polarization. (Some confirmation of the prediction is observed with 
Cdc42 manipulation in yeast, where the frequency of spontaneous polarization is inversely dependent on 
the amount of Cdc42 [T].) Replacing a cytosolic protein by a fusion protein with lower mobility has been 
experimentally done in budding yeast [5] . Our results show that reducing D for a fixed e may lead to loss of 
polarity, as the border between wave-pinning and homogeneous regimes is shifted (e.g. sec Fig. 14.21 second 
panel). Furthermore, Rho protein cycling between membrane and cytosol is affected by proteins called 
GDIs. We have previously shown that the time spent in the cytosol vs membrane affects the effective 
diffusion coefficient of the inactive form D v [T21 [T5] , which thus affects the dimensionless parameter D. 
This suggest that regulation of the GDIs is yet another possible mechanism for regulating polarity [3J). 
Experiments in budding yeast show that knock down (i.e., replacement with a non- functional version) of 
GDI coupled with treatment to disable a second redundant Cdc42 membrane recycling pathway leads to 
rapid dissipation of polarity 33 . This supports our predictions. 
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6. Appendix. 

6.1. Integral Reduction at Finite D. We perform an integral reduction of (|4.1[) for finite D, 
similarly to the treatment of the D — >• oo case in Section 14.21 We view (|4.1|) as an ODE with x as the 
"time" variable. First, add (|4.1a .b) to obtain 

e 2 u xx + Dv xx = 0. (6.1) 

Integrate this equation twice and use the no-flux boundary conditions to obtain 

e 2 

— u + v = A, (6.2) 

where A is an integration constant. Solving the above for v and substituting this into (|4.1a|) . we reduce 
the system to a single equation for u: 

= e 2 u xx + f D (u,A), (6.3) 



where fu(u,A) = f(u,A— e 2 u/D). The rest follows along exactly the same lines as in Section l4~2l 

We note two differences. Recall that the function /(it, v) is bistable in u for fixed v satisfying 
"min < v < "Umax- We can see from (|6.2[) that if 

V m in < A < (6.4) 

then fri(u,A) is bistable in u (for u in a finite range) for e small enough, assuming that f(u,v) is a 
smooth function of v. In this case, the function: 



F D (u,A,B) = -B- f f t 
Jo 



(s,A)ds (6.5) 



will have the form of a double well potential, whose local minima correspond to the stable zeros of the 
bistable function /d(u, A). This restriction on the size of e was absent in the D — > oo case. This can also 
be seen by formally taking the limit as D — > oo in (|6.2[) . which yields A = v. 
The integral conditions (|4.8p and (|4.9|) . in the finite D case have the form: 

l = ef - dU , (6.6) 

„ / e 2 \ f" 1 udu . . 

if = A+e 1-- - 6.7 

V DjJ U0 y/F D (u,A,By 

where uo < u% are the two middle roots of the equation Fd(u, A, B) — 0. It is easy to see that these 
conditions reduce to (|4.8I) and (|4.9[) in the D —> oo limit. One difficulty here is that it is not possible 
to eliminate e to obtain a relation analogous to (|4.10j) . since Fd(u,A,B) has an e dependence. It is 
nonetheless possible to use the above as a basis for a continuation algorithm, and we have seen that the 
results using these relations match with those obtained using a direct discretization of (|4.ip [TT] . 

6.2. Derivation of the Expression for K r . In this appendix, we derive expressions (|4.25[) and 
(14.261) . Consider (|4.4[) when (|2.17|) is used for the reaction term where a > 1. Take any < K < K q 
where K q is given in (|4.23[) . Let v$ be the middle root of equation (|4.22[) and let uq = u m (vo) (note 
that we referred to vq as and uq as in Section |4.3[) . We saw that (uq,vo) is a stable spatially 
homogeneous solution to (|4.3[) for 

e > e = ; 1 (6.8) 

^(l + Ho) 2 ) 

where we used (|4.24p . At e = eo, we demonstrated computationally that we have a pitchfork bifurcation. 
We now perform a perturbation calculation to study this bifurcation (see, for example, [13] or [S]). 

Let us restate our problem for future reference. For algebraic convenience, we shall work with A = 1/e 2 
instead of e. We study the bifurcation of the steady state solution (u, v) = (u ,v ) of the system 

— -\(u 2 -l)[u+^===) =0, (6.9) 



9x 2 y yf\ + (av) 

du 
dx 



udx = K, — 



0, (6.10) 



x=0,l 
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at the bifurcation point A = Aq. The values Aq and uq can be expressed in terms of vq: 



where vq is the middle root of: 



A = ^ 2 (1 + (av ) 2 ), u Q = 7= ™° = , (6.11) 

Vi + {av y 



vo - ,7 . a = K. (6.12) 
^1 + (av ) 2 



Note that i>o can thus be viewed as a function of K where < K < K q . It is easy to see that vq(K) is a 
decreasing function of K. As K ranges from to K qi vq ranges from to — v a 2 / 3 — 1/a. 
We introduce a small parameter <5 and seek a solution of the form: 

u = u Q + 5ui + S 2 u 2 + S 3 u 3 H . (6.13) 

We introduce a similar expansion for A and v. Substitute these into (|6.9p and (|6.10[) and collect like 
terms in S. The 0(1) relation is just Xo,uq,vq substituted into (|6.9|) and (|6.10|) . and thus does not give 
us anything interesting. At 0(<5) we have: 

d 2 { A f ^ 

^ " Mu 2 - 1) + (1 + (fl ° o)2)3/2 ^) = 0,-i + / tud* = 0. (6.14) 
For a function / defined on < x < 1, define: 

Qf = f-A ( fdx, A= a (6.15) 

Jo (i + (aw rr /2 

Using this and (|6.11[) . we may rewrite (|6.14j) as: 

d 2 ui 



dx 2 



7r 2 Q Ul =0, (6.16) 



where we have Neumann boundary conditions at x = 0, 1. The only nontrivial solutions to the above are 
constant multiples of cos(7rx). We thus let: 

U\ = cos(-kx). (6-17) 

Other choices of u\ merely amounts to a rescaling of 8. Note that v\ = by (|6.14|) . 
At 0(<5 2 ) we have, after some simplification: 

2 u 

+ n 2 Qu 2 = 2A u u 2 + Ai(ug - l)m. (6.18) 

Given that the operator on the left hand side is self-adjoint with the null space spanned by cos(7rx), we 
require that the left hand side be orthogonal to this. From this, we easily conclude: 

Ai=0. (6.19) 

We may then solve for u 2 imposing orthogonality with respect to u\ = cos(7rx) to obtain: 



u 2 =a + l3cos(7rx), a= — ^ ± , /? = ^av y/l + (av ) 2 , (6.20) 

where A was given in (|6.15p . 
At 0(<5 3 ) we have: 

+ 7I ' 2 2 U 3 = A 2 (uq - l)«i + 2X u uiQu2 + X (2uqu 2 + u\)ui, (6-21) 

where we have used Ai = 0. The left hand side must be orthogonal to ui, from which we obtain the 
following expression for A 2 : 

A 2 = A f/ (l- u 2 )u 2 dx) ( f (2u u 2 Qu 2 + (2u u 2 + u 2 )u 2 )dx] . (6.22) 
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Given that A = Ao + 8 2 \-2 + • • • , the sign of A2 determines whether the pitchfork bifurcation is subcritical 
or supercritical. We thus seek the point at which A2 changes sign as K is varied. The sign of A2 is 
determined by the sign of the last integral in (|6.22[) . This integral can be computed as: 

f 1 3 
1=1 (2u u 2 1 Qu 2 + {2u u 2 + u 2 1 )u 2 1 )dx = u a(2- A) +u /3 + -. (6.23) 
Jo 8 

We may simplify this expression to find: 

I = (s 2 - 1) (| + ^) + I * = Vl + W 2 - (6.24) 

Using properties of vq (K) , we see that s is an increasing function of K and varies between 1 < s < ^/a 
for < K < K C] . Dividing the above by s 2 — 1, we obtain the left hand side of (|4.25[) . which is monotone 
in s for 1 < s < -^fa. We see that there is only one value of s and hence K at which / changes sign. This 
is the value of K r we seek. 
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